############################
############################
##### Replication code #####
##### for ##################
##### Still 'the Domain ####
##### of Men'? #############
##### Weeks and Masala #####
############################
############################

rm(list = ls())

### set working directory 

## load the bill level data 

load("wm_bill_level.RData")
head(dat)
nrow(dat) ##2366 

## load the region-year level data 

load("wm_region_level.RData")
head(test)
nrow(test) #22

### load the individual councillor level data 
load("wm_individual_level.RData")
head(billdat)
nrow(billdat) # 238

options(scipen=999)

#################################################################
##################### Install packages if necessary #############
#################################################################

#install.packages("ggplot2")
#install.packages("gmodels")
#install.packages("stargazer")
#install.packages("aod")

library(ggplot2)
library(gmodels)
library(stargazer)
library(aod)


#################################################################
#################################################################

## How many bills come from Councillors? 
table(dat$P_C, dat$region)
table(dat$P_C[dat$region=="Cam"]) ##72% of bills proposed by councillor
CrossTable(dat$P_C[dat$region=="Cam"], prop.t=TRUE, prop.r=TRUE, prop.c=TRUE)

table(dat$P_C[dat$region=="Cal"]) ##71% of bills proposed by councillor
CrossTable(dat$P_C[dat$region=="Cal"], prop.t=TRUE, prop.r=TRUE, prop.c=TRUE)

## 1667 bills come from councillors, of these 303 from women 
table(dat$P_C, dat$P_W)

table(dat$P_C[dat$region=="Cal"], dat$P_W[dat$region=="Cal"]) ## In Calabria only 27 bills proposed by women (3.5%)
table(dat$P_C[dat$region=="Cam"], dat$P_W[dat$region=="Cam"]) ## In Campania 276  bills proposed by women (30%)

### How many of the bills that pass come from councillors ? 
table(dat$passed)
table(dat$passed, dat$P_C)

table(dat$passed[dat$region=="Cal"], dat$P_C[dat$region=="Cal"]) ##57% in Calabria
CrossTable(dat$passed[dat$region=="Cal"], dat$P_C[dat$region=="Cal"], prop.t=TRUE, prop.r=TRUE, prop.c=TRUE)

table(dat$passed[dat$region=="Cam"], dat$P_C[dat$region=="Cam"]) ##44% in Campania 
CrossTable(dat$passed[dat$region=="Cam"], dat$P_C[dat$region=="Cam"], prop.t=TRUE, prop.r=TRUE, prop.c=TRUE)


#################################################################
################### Replicates Figure 1 ########################


p2 <- ggplot(test, aes(year, share_women_proposed, group = Region, colour = Region))  + geom_point(size = 2) + geom_line(size = 0.8) +
     theme(legend.position = "none") + 
theme_bw() +
theme(text=element_text(color="black"))   + xlab("") + ylab("Share bills sponsored by women")

p2<- p2 + scale_x_continuous(breaks=seq(2007,2017,1), limits=c(2007,2017)) + scale_y_continuous(limits=c(0,65)) + theme_classic() + 
    theme(axis.title=element_text(size=15)) + geom_vline(xintercept=c(2010), linetype="dashed")  + scale_colour_grey(start=0.8, end=0.2)


p2 + geom_line(linetype="dotted", size = 0.8, aes(x=year,y=mean_womenc, group = Region, colour = Region)) + 
  annotate("text", x=2016, y=23, label= "% Women in Council, Campania") + 
  annotate("text", x = 2016, y=3, label = "% Women in Council, Calabria")

## ggsave("percent_bills_proposed_by_women.tiff", width = 12, height = 7)


#############################################################
#############################################################
####################
### Table 1: DV share of bills proposed 
####################

table(test$share_women_proposed[test$region=="Campania"], test$year[test$region=="Campania"])
table(test$share_women_proposed[test$region=="Calabria"], test$year[test$region=="Calabria"])

table(test$mean_womenc[test$region=="Campania"], test$year[test$region=="Campania"])

###DV: share of all bills proposed by councilors with a woman sponsor or co-sponsor. 

mod3a<-lm(share_women_proposed ~ quota   + as.factor(region) + as.factor(year), data = test)
summary(mod3a)

stargazer(mod3a,
          font.size = "small",
          title="Regression Results",
          align=TRUE, dep.var.labels=c(
            "% Bills Proposed by Women"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="Table1a.doc")

#######################################################
###### Individual councillor level data ###############
#### Mean number of bills proposed per year pre vs post quota by men and women in Campania: 
### For this analysis, use councillors elected 2005, 2010, 2015
#######################################################

head(billdat)

billdat$quota<-0
billdat$quota[billdat$COUNCIL.TERM>2005]<-1
table(billdat$quota)

billdat$prepost<-NA
billdat$prepost[billdat$COUNCIL.TERM==2005]<-"pre"
billdat$prepost[billdat$COUNCIL.TERM==2010]<-"post1term"
billdat$prepost[billdat$COUNCIL.TERM==2015]<-"post2term"
table(billdat$prepost)

billdat$woman<-0
billdat$woman[billdat$GENDER..W.1.==1]<-1
table(billdat$woman)

billdat2<-billdat[ which(billdat$COUNCIL.TERM<2016), ]

mean(billdat2$BILLS_PER_YEAR[billdat2$quota==0])
mean(billdat2$BILLS_PER_YEAR[billdat2$quota==1])

## mean pre to post quota : women // 1.94 increase
t.test(billdat2$BILLS_PER_YEAR[billdat2$quota==0 & billdat2$woman==1], billdat2$BILLS_PER_YEAR[billdat2$quota==1 & billdat2$woman==1] )

### and a one-sided test 
t.test(billdat2$BILLS_PER_YEAR[billdat2$quota==0 & billdat2$woman==1], billdat2$BILLS_PER_YEAR[billdat2$quota==1 & billdat2$woman==1],
   var.equal = TRUE, alternative = "less" )

## men // 1.151 increase
t.test(billdat2$BILLS_PER_YEAR[billdat2$quota==0 & billdat2$woman==0], billdat2$BILLS_PER_YEAR[billdat2$quota==1 & billdat2$woman==0],
  var.equal = TRUE, alternative = "less" )

## mean differences between men and women from 2005 
t.test(billdat2$BILLS_PER_YEAR[billdat2$prepost=="pre" & billdat2$woman==0], billdat2$BILLS_PER_YEAR[billdat2$prepost=="pre" & billdat2$woman==1],
   var.equal = TRUE, alternative = "less"  )
## from 2010
t.test(billdat2$BILLS_PER_YEAR[billdat2$prepost=="post1term" & billdat2$woman==0], billdat2$BILLS_PER_YEAR[billdat2$prepost=="post1term" & billdat2$woman==1],
   var.equal = TRUE, alternative = "less"  )
## from 2015
t.test(billdat2$BILLS_PER_YEAR[billdat2$prepost=="post2term" & billdat2$woman==0], billdat2$BILLS_PER_YEAR[billdat2$prepost=="post2term" & billdat2$woman==1],
   var.equal = TRUE, alternative = "less"  )


################################################
######## Table 2: gender quota and proposals by governing party status 
#################################################

### look at pre-post bill proposals from women across different types (govt, opp, bipartisan, both)
##use data only from campania 

### How many bills in each category? 
table(dat$Bill_Proposer_Bipartisan)
table(dat$Bill_Proposer_Opposition)
table(dat$Bill_Proposer_Govt)

camdata<-test[ which(test$region=='Campania'), ]

mod3ab<-lm(share_proposed_women_govt ~ quota, data = camdata)
summary(mod3ab)
mod3ab2<-lm(share_proposed_women_opp ~ quota  , data = camdata)
summary(mod3ab2)
mod3ab3<-lm(share_proposed_women_bipartisan ~ quota  , data = camdata)
summary(mod3ab3)

stargazer(mod3ab,mod3ab2, mod3ab3,
          font.size = "small",
          title="Regression Results",
          align=TRUE, dep.var.labels=c(
            "% Bills Proposed by Women, Governing parties",
            "% Bills Proposed by Women, Opposition parties", 
            "% Bills Proposed by Women, Bipartisan"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="Table2_parties.doc")


############################
############################
########## Table 3: Effectiveness of women councillors (Bills passed)
############################
############################


dat2<-dat[dat$P_C==1,]

### DV = binary variable indicating whether the bill proposed was passed or not. Sample includes only those bills proposed by councillors.
## N = 1688

nrow(dat2)

mylogit <- glm(passed ~ quota + P_W + quota*P_W + as.factor(region) + as.factor(year), data = dat2, family = "binomial")
summary(mylogit)

OR.vector <- exp(mylogit$coef)
CI.vector <- exp(confint(mylogit))
p.values <- summary(mylogit)$coefficients[, 4]

stargazer(mylogit, coef = list(OR.vector), ci = T, 
          ci.custom = list(CI.vector), p = list(p.values), 
          single.row = T, type = "html", out="Table2.doc")

##################
### Logit regression / Odds ratio models of bill passage by governing party status
##################

#### Look at govt bills only 
dat3<-dat[ which(dat$Bill_Proposer_Govt==1), ]
nrow(dat3)
table(dat3$region) ##465 bills proposed by governing party members in campania 

mylogit2 <- glm(passed ~ quota + P_W + quota*P_W , data = dat3, family = "binomial")
summary(mylogit2)

OR.vector2 <- exp(mylogit2$coef)
CI.vector2 <- exp(confint(mylogit2))
p.values2 <- summary(mylogit2)$coefficients[, 4]


### same for opposition 

dat4<-dat[ which(dat$Bill_Proposer_Opposition==1), ]
nrow(dat4)
table(dat4$region) ##465 bills proposed by governing party members in campania 

mylogit3 <- glm(passed ~ quota + P_W + quota*P_W , data = dat4, family = "binomial")
summary(mylogit3)

OR.vector3 <- exp(mylogit3$coef)
CI.vector3 <- exp(confint(mylogit3))
p.values3 <- summary(mylogit3)$coefficients[, 4]


stargazer(mylogit2, mylogit3, coef = list(OR.vector2, OR.vector3), ci = T, 
          ci.custom = list(CI.vector2, CI.vector3), p = list(p.values2, p.values3), 
          single.row = T, type = "html", out="Table_Appendix_odds_ratios.doc")

########################################################
########################################################
### Table 4 ############################################
### Leadership roles ###################################
########################################################

#### test for significant differences in women's leadership before and after quota law

coun <- billdat
##read.csv("council_dat.csv")

coun$reelect<-0
coun$reelect[coun$RE.ELECTION..FOR.A.FOLLOWING.TERM.==1]<-1
table(coun$reelect)

coun$woman<-0
coun$woman[coun$GENDER..W.1.==1]<-1
table(coun$woman)

coun$quota<-0
coun$quota[coun$COUNCIL.TERM>2005]<-1
table(coun$quota)

coun$leadpc<-0
coun$leadpc[coun$LEADERSHIP.PERMANENT.COMMITTEE==1]<-1
table(coun$leadpc)

coun$leadsc<-0
coun$leadsc[coun$LEADERSHIP.SPECIAL.COMMITTEE==1]<-1
table(coun$leadsc)

coun$vomm<-0
coun$vomm[coun$V.COMMITTEE.MEMBERSHIP==1]<-1
table(coun$vomm)

coun$speakersoffice<-0
coun$speakersoffice[coun$SPEAKER.S.OFFICE==1]<-1
table(coun$speakersoffice)

coun$year<-coun$COUNCIL.TERM

### lead perm committee 

coun3<-coun[which(coun$leadpc==1), ]
nrow(coun3)

coun3$preq<-0
coun3$preq[coun3$year==2005]<-1
table(coun3$preq)

mean(coun3$woman[coun3$preq==1])
sd(coun3$woman[coun3$preq==1])

mean(coun3$woman[coun3$preq==0])
sd(coun3$woman[coun3$preq==0])

chisq.test(table(coun3$preq, coun3$woman))

### lead special committee 

coun3<-coun[which(coun$leadsc==1), ]
nrow(coun3)

coun3$preq<-0
coun3$preq[coun3$year==2005]<-1
table(coun3$preq)

mean(coun3$woman[coun3$preq==1])
sd(coun3$woman[coun3$preq==1])

mean(coun3$woman[coun3$preq==0])
sd(coun3$woman[coun3$preq==0])

chisq.test(table(coun3$preq, coun3$woman))

#### 5th committee member

coun3<-coun[which(coun$vomm==1), ]
nrow(coun3)

coun3$preq<-0
coun3$preq[coun3$year==2005]<-1
table(coun3$preq)

mean(coun3$woman[coun3$preq==1])
sd(coun3$woman[coun3$preq==1])

mean(coun3$woman[coun3$preq==0])
sd(coun3$woman[coun3$preq==0])

chisq.test(table(coun3$preq, coun3$woman))

### office of speaker 

coun3<-coun[which(coun$speakersoffice==1), ]
nrow(coun3)

coun3$preq<-0
coun3$preq[coun3$year==2005]<-1
table(coun3$preq)

mean(coun3$woman[coun3$preq==1])
sd(coun3$woman[coun3$preq==1])

mean(coun3$woman[coun3$preq==0])
sd(coun3$woman[coun3$preq==0])

chisq.test(table(coun3$preq, coun3$woman))
fisher.test(table(coun3$preq, coun3$woman))


#############################################################
#################### Replicates Figure 2 ####################

p3 <- ggplot(test, aes(year, share_genderpol_proposed, group = Region, colour = Region))  + geom_point(size = 2) + geom_line(size = 0.8) +
     theme(legend.position = "none") + 
theme_bw() +
theme(text=element_text(color="black"))   + xlab("") + ylab("Share of bills related to women's interests")

p3<- p3 + scale_x_continuous(breaks=seq(2007,2017,1), limits=c(2007,2017)) + scale_y_continuous(limits=c(0,20)) + theme_classic() + 
    theme(axis.title=element_text(size=15)) + geom_vline(xintercept=c(2010), linetype="dashed")  + scale_colour_grey(start=0.8, end=0.2)


p3
## ggsave("percent_gender_bills.tiff", width = 12, height = 7)


##############################
####### Table 5: Substantive Rep ##
##############################

## DV is the share of bills proposed out of all bills proposed by councillors

mod3<-lm(share_genderpol_proposed ~ quota + funemp  + as.factor(region) + as.factor(year), data = test)
summary(mod3)## quotas lead to a significant increase in the s

### confirm that results hold without women's unemployment
mod3a<-lm(share_genderpol_proposed ~ quota   + as.factor(region) + as.factor(year), data = test)
summary(mod3a)## quotas do not lead to a significant increase 

test$quota_adopt<-0
test$quota_adopt[test$region=="Campania" & test$year>2008]<-1
table(test$quota_adopt, test$year)

### confirm that results are the same when considering quota adoption rather than implementation
mod3a<-lm(share_genderpol_proposed ~ quota_adopt   + as.factor(region) + as.factor(year), data = test)
summary(mod3a)## quotas do not lead to a significant increase 

stargazer(d3,
          font.size = "small",
          title="Table 4: Effect of Quota Law on Substantive Representation of Women�s Priorities",
          align=TRUE, dep.var.labels=c(
                        "Share Women's Issues Bills"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="Table4a.doc")


###########################
##### Appendix
###########################

### Figure A1

p <- ggplot(test, aes(year, mean_womenc, group = Region, colour = Region))  + geom_point(size = 2) + geom_line(size = 0.8) +
     theme(legend.position = "none") + 
theme_bw() +
theme(text=element_text(color="black"))   + xlab("") + ylab("Share women in regional council")

p<- p + scale_x_continuous(breaks=seq(2007,2017,1), limits=c(2007,2017)) + scale_y_continuous(limits=c(0,50)) + theme_classic() + 
    theme(axis.title=element_text(size=15)) + geom_vline(xintercept=c(2010), linetype="dashed")  + scale_colour_grey(start=0.8, end=0.2)

p

## ggsave("women_on_council.pdf", width = 12, height = 7)

##############################################
##### Table A1 ###############################
### Effect of quota law on women in council 

mod3b<-lm(mean_womenc ~ quota   + as.factor(region) + as.factor(year), data = test)
summary(mod3b)## quotas lead to a significant increase in the s


stargazer(mod3b,
          font.size = "small",
          title="Regression Results",
          align=TRUE, dep.var.labels=c(
            "% Women in Council"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="TableA1.doc")

############ Table A2: Drop years 2016 and 2017

mod3ab<-lm(share_women_proposed ~ quota   + as.factor(region) + as.factor(year), data = test, year < 2016)
summary(mod3ab)

stargazer(mod3ab,
          font.size = "small",
          title="Regression Results",
          align=TRUE, dep.var.labels=c(
            "% Bills Proposed by Women"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="TableA2.doc")

#############################
########### Figure A2
### Histograms pre and post quota 
#############################

hist(billdat$BILLS_PER_YEAR[billdat$quota==0 & billdat$woman==1], main = "Pre-quota law",
     xlab = "Average number of bils sponsored per year by women councillors",  ylim = c(0,10))
abline(v = mean(billdat$BILLS_PER_YEAR[billdat$quota==0 & billdat$woman==1]),                       # Add line for mean
       col = "red",
       lwd = 1)
text(x = 2.7,                   # Add text for mean
     y = 8,
     paste("Mean =", 2.1),
     col = "red",
     cex = 1)

hist(billdat$BILLS_PER_YEAR[billdat$quota==1 & billdat$woman==1], main = "Post-quota law",
     xlab = "Average number of bils sponsored per year by women councillors",  ylim = c(0,10))
abline(v = 4.1,                       # Add line for mean
       col = "red",
       lwd = 1)
text(x = 5.2,                   # Add text for mean
     y = 8,
     paste("Mean =", 4.1),
     col = "red",
     cex = 1)

#################################
#### Table A3: 
#### quota effects from adoption
#################################  

test$Campania<-0
test$Campania[test$region=="Campania"]<-1
table(test$Campania) 
table(test$region)

test$yearf<-as.factor(test$year)

test$quota_adopt<-0
test$quota_adopt[test$region=="Campania" & test$year>2008]<-1
table(test$quota_adopt, test$year)

mod3d<-lm(share_women_proposed ~ quota_adopt   + as.factor(region) + as.factor(year), data = test)
summary(mod3d)

stargazer(mod3d,
          font.size = "small",
          title="Regression Results",
          align=TRUE, dep.var.labels=c(
            "% Bills Proposed by Women"),
                                       omit.stat=c("LL","ser","f"), 
          no.space=TRUE,
          omit = c("year", "region"),
          column.sep.width = "-15pt", type="html", out="TableA3_Appendix.doc")












